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Abstract 


During geometry processing , tangent directions at the data points are 
frequently readily available from the computation process that generates 
the points. It is desirable to utilize this information to improve the 
accuracy of curve fitting and to improve data reduction. 

This paper presents a curve fitting method which utilizes both posi- 
tion and tangent direction data. This method produces G 1 nonrational 
B-spline curves. From the examples , the method demonstrates very good 
data reduction rates while maintaining high accuracy in both position 
and tangent direction. 


1 Introduction 


Path tracing is a common technique utilized by many numerical processes 
in computers, e.g., computing the intersection between surfaces, evaluating 
the silhouette curves, and solving partial differential equations. Oversampling 
is a method frequently used by tracing algorithms to overcome uncertainty in 
tracing step selection and to ease various stability problems. A large amount of 

*This work was performed for NASA under Contract NAS 2-12961 
iTo be submitted as a NASA Contractor Report 
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points results from tracing, partially due to oversampling. For the convenience 
of subsequent data handling and storage, it is desirable to fit the data by curves 
with less data. Hence data reduction is an important aspect of fitting traced 
data. In this paper, the word fitting means to pass a curve close to the data 
but not necessarily exactly through the data. 

Tangent direction at the data points is usually used to determine the marching 
direction in path tracing; hence it is the by-product of the tracing algorithms. 
Traditional curve fitting algorithms seldom take advantage of this tangent 
information to improve the accuracy of curve fitting or to enhance data reduc- 
tion. 

Robust tracing algorithms need to detect discontinuity in the path curve. 
Hence, the resulting points can be grouped into pieces forming continuous 
paths. Thus, in this paper we assume the given data set forms a continuous 
path. 

The problem at hand is different than fitting inaccurate data, that is data with 
relatively large error, e.g., digitized data. In such cases, the fitting curve is 
sought to best approximate the shape, not the individual points. High order 
information such as tangent is meaningless in such case. However, in tracing, 
the data, both position and tangent direction, are usually considered to be 
“exact,” and the fitting curve is required to be as close to the data as possible. 


Cubic spline fits or least-squares spline fits have been traditionally used for 
fitting B-spline curves [1, 2, 3, 4, 5, 10]. Cubic spline fits always produce 
more data than the input. To achieve a given tolerance with the least-squares 
methods, a fit-then-test procedure with an increasing number of control points 
is necessary. When the number of control points is increased to the number of 
input data, least-squares methods converge to interpolation. According to the 
author’s experience, least-squares methods produce reasonable results when 
the required fitting accuracy is low (> 10 -2 ). For slightly higher accuracy 
(< 10 2 )> least-squares methods always result in interpolation, producing no 
data reduction. For fitting with B-splines, interpolation produces a curve with 
more data than the original if the knots are counted. Various data reduction 
schemes [11, 12, 13] have been suggested by other authors. However it is very 
difficult to make any reduction with these methods when the required accuracy 
is high. 

Kallay [7] presented a method to fit a G 1 polynomial curve to a planar set of 
points with tangent directions. Piegl [6] and Chou [8, 9] also devised methods 
to fit C 1 rational B-spline curves to the same data for planar and for three- 
space cases. However, none of these methods have achieved high data reduction 
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rates. Besides, there are systems which are based on integral B-splines, hence, 
can not utilize the methods in [8, 9]. 

In this paper, we describe a method which follows the same divide-and-conquer 
strategy used in [6, 8, 9] but produces an integral B-spline fitting curve. A 
new tolerance checking method is presented to speed up the tolerance checking 
process. The data reduction rate is compared favorably with that of previous 
methods. 

The method is as follows: first, find the longest run of data that can be fit with 
one cubic Bezier curve, starting from the last fit data point. At the beginning, 
the longest run starts from the first data point. This process repeats until all 
the data are fit. The Bezier curves from the fits axe then connected together 
to form a G 1 B-spline curve. 

To locate the longest run of data that can be fit by one cubic Bezier curve, 
an attempt is made to fit the given list of data. If the data can not be fit by 
one Bezier curve with the given tolerance, the set of data is reduced, and the 
process to fit one cubic Bezier curve is repeated recursively. If the data can 
be fit by one Bezier curve, the set of data is enlarged, and the process to fit 
one cubic Bezier curve is repeated recursively. The recursion stops when the 
the longest run of data that can be fit is found. 

Given a list of data, the following steps are used to obtain the Bezier fitting 
curve: 


1. The first and last control points of the Bezier curve are the first and last 
points of the data, respectively. 

2. The tangents at the first and last points of the Bezier curve are in the 
tangent directions of the first and last points of the data, respectively. 

3. A cubic Bezier interpolation curve is created for each of the inner data 
points (data except the first and last point). 

4. A cubic Bezier fitting curve is created by weighting the interpolation 
curves. 

5. Checks are performed to insure the fitting curve is within the given tol- 
erance. 


In Step 3, the interpolation curve has the same control points and tangent di- 
rections as the fitting curve. In addition, the interpolation curve goes through 
the inner data point. A planar case happens when all of the first and last 
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points, the end tangent directions, the inner data point, and the tangent di- 
rection at the inner data point are on a plane. In this case, the tangent 
direction of the interpolation curve at the inner data point also agrees to the 
given tangent direction at the point. 

In the following sections, Steps 3, 4 and 5 are discussed in detail. 


2 Interpolation Cubic Computation 


In this section we discussed the methods to obtain the interpolation curves. 
Most of the equations in this section can also be found in [8, 9]. 

A cubic Bezier curve can be written as: 

C{u) = P 0 {BZ{u) + Bf(u)) + atoB^u) + ^B^u) + P 3 (B*(u) + £f(tz)) (1) 

where 


. B*(u) = (JW(1 - u) 3 - ' are the cubic Bernstein polynomials; 

• to and t\ are the end tangent directions (unit length); and 

• Po, Psi P\ — Pq + Cato., P 2 = P 3 + are the control points. 

The above equation assumes, without loss of generality, that the cubic starts 
at u = 0 and ends at u = 1. An illustration of the variables is in Figure 1. 

For a given data point P, the interpolation curve goes through this point. 

P = P 0 (Bq(u) + B*(u)) + atoBfiu) + /?t x £ 2 3 (u) + P 3 (P 2 3 («) + P|(«)) (2) 

There are three variables in Equation (2): a, j3, and u. Equation (2) actually 
represents three equations, for the z, y , and z components, respectively. When 
to, tj, and P0P3 span the three-space, it is proved in [8] that there exists at 
most one set of solution (a, /?, u ) with u € [0,1]- u can be found by performing 
dot products on Equation (2) with < 3 = f 0 x tj: 

= (3) 
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This is a cubic Equation in u. Exactly one solution exists in [0, 1], if and only 
if € [0> !]• Once u is obtained, a and /?, and hence the interpolation 

curve, can be found from Equation (2). 

When to, ti, and PqPo are coplanar, Equation (2) has no solution if P is not 
in the plane. If P is in the plane, Equation (2) provides only two independent 
equations. If the tangent direction at P, denoted as <2 (unit length), is also 
in the plane, an additional equation is provided by forcing the interpolation 
curve to agree to ti at P. 

«2 = (ft - ft)B, 2 (i) + ct„(Bl(u) - B\(u)) + pU(B\(u) - B|(S)), (4) 

where 8 is a multiplier, and Bf(u) = ^u*(l— u) 2_l are the quadratic Bernstein 
polynomials. 

Equations (2) and (4) provide four equations and have unknowns: (a, /?, 8, u). 
The equations can be combined and simplified to form a cubic equation in 
u [9]. 

- [((P - P 0 ) • h)(U ■ to) + ((P - Po) • to)(ti ■ h)] - (5) 

3[((P — Po) • ti)(to • ?i)]i2 + 

3[((p3 — Po) • ti)(to ■ ii)\v? + 

[((P 3 — Po) • ^l )(*2 ' to) + ((P 3 — Po) • t 0 )(t 2 ■ ^l)]u 3 = 0, 

where 

Jc is the normal of the plane, <0 = k x t 0 , h = k x t\, t 2 = k x t 2 . (6) 


Once u is obtained, a, /3, and 5 can be computed. From a and /? the in- 
terpolation Bezier curve can be obtained. However, in the planar case, there 
may be more than one solution set (a, (3, 6, u) with u G [0, 1). Even with the 
restrictions that the tangent directions of the Bezier curve should agree to the 
given data, i.e., a > 0,/? < 0, and 8 > 0, multiple solutions may still occur. In 
addition, Equations (2) and (4) do not have solutions in some special cases, 
for example, when to,t\, and t 2 are parallel. 


3 Fitting Curve Determination 


When multiple solutions occur in computing the interpolation curve, the so- 
lution that has the closest matching end tangents with those from the neigh- 
boring data point is chosen as the interpolation curve. On the other hand, a 
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valid interpolation curve may not exist. In such a case, the data point does 
not participate in determining the fitting curve (Step 4). 


The fitting curve is obtained by weighting the interpolation curves. Since all 
the interpolation curves have the same start and end points and the same start 
and end tangent directions, locating the two inner control points (Pi and P 2 ) 
is enough to fix the fitting curve. This is equivalent to determining the a and 
/? of the fitting curve, denoted as a and j3. We express d as a weighted sum 
of as from the interpolation curves, denoted as a,. 


a = 


Ei w.-Qi 
E. m 


(7) 


Three possible methods to compute d are listed below. can be computed by 
similar methods. 


1 : Wi = 1 . 

2: wj = 1, when aj = max or min a, ; w t = 0, for the rest of i. 

3: Wi = Pf(u,), where ti, is the parameter of the interpolation curve at the 
ith data point. 


Method 1 gives even weights to all the interpolation curves. Roughly speak- 
ing, Method 2 takes the average of the maximum and minimum interpolation 
curves. The third method gives heavier weight to the interpolation curve whose 
data point is closer to P 0 . This results in a fitting curve that deviates less from 
the data point when the point is closer to Pq. Since P\ has greater influence 
than P 2 on the part of the curve closer to Pq , Method 3 reduces the error of 
the fitting curve. Therefore, Method 3 is used in this paper. The formula 
Wi = B$(iii) is used to compute /3. 


4 Tolerance Handling 


Since we are interested in achieving high accuracy, we would like to check the 
distance between the data point and the fitting curve. However, it is quite 
expensive to compute the closest distance from a point to a curve. In order 
to reduce this cost, we first perform two pre- tests: one gives a quick rejection 
of the fitting curve, the other raises a quick acceptance of the fit. If the fit 
falls through both of the tests, the more costly closest distance computation 
is performed. 
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When we compute the interpolation curves for a list of data, we keep the 
maximum and minimum values of the u>,a, (and w t 0,) for the curves obtained 
so far. When the difference between the maximum and minimum is greater 
than a constant value, we declare the fit to be unsuccessful. This condition 
may be reached before all the data points are processed, and usually it happens 
during the first several recursive runs of fitting. Hence, this test serves as a 
quick rejection to data sets that variate too much to be fit by a cubic, and 
it significantly accelerates the process in locating the longest run of data. 
A good value for the constant that determines the rejection is two order of 
magnitude of the given tolerance for a tolerance of 10 -4 . This formula is purely 
experiential. Note that this method may declare a fit to be unsuccessful, even 
though the fitting curve is within the given tolerance. 

For the quick acceptance test, we observe that 

\C-P\ < \C{u)-P\ = lAaB^uj + A^u)! < |Aa|fi 1 3 (u) + |A/3|H 2 3 (t2), (8) 

where Aa — a — a and A 0 = 0-/3. Both i? 3 (u) and B * («) are available from 
the process computing the interpolation curve. Therefore, for each point, this 
estimation of the upper bound for the error can be obtained by performing 
two multiplications, two subtractions and one addition. 

Finally, computing the closest distance between the fitting curve and a data 
point is necessary only if the data point fails the quick acceptance test. When 
this happens, we use the Newton-Raphson method to locate the closest point 
on the curve to the data point. The parameter u of the interpolation curve 
at the data point serves as a very good initial guessing point for the Newton- 
Raphson method. 


5 Examples 


In this section, we show the results from applying the above fitting method 
to four data sets. The first data set (Half-Circle) is created by sampling 101 
points on a half circle. The second data set (Torus-Plane) with 361 points 
is created by tracing the intersection curve of a torus and a plane. Both of 
these data sets are planar. But only the first data set is in the x-y plane. The 
third data set (Planar-Spiral) contains a planar section and a spiral three-space 
section with a total of 51 data points. The remaining data set (Mixed-Spiral) 
has 101 data points with mixed planar and three-space sections. Figures from 
the sample fits are shown in Figures 2 to 5. The first three figures can be 
compared with those from the rational curve fits in [9]. The last figure can be 
compared with that in [8]. 
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Half- 

Circle 

Torus- 

Plane 

Planar- 

Spiral 

Mixed- 

Spiral 

Tolerance 

10" 3 

10" 6 

io - 2 

10" 4 

2 • 10" 1 

10' 3 

10" 5 

10“ 3 

10“ 5 

From 

(101) 

(361) 


(51) 


(101) 

Input 

303 

1805 


255 


505 

Current 

(10) 

(22) 

(34) 

(61) 

(13) 

(34) 


(64) 

(85) 

Method 

34 

70 

140 

248 

56 


164 

260 

344 

Method 

(10) 

(58) 

(43) 

(112) 

(16) 

(37) 

(61) 

(82) 

(100) 

in [9] 

44 

236 

219 

564 

84 

189 


414 

504 

Inter- 

(152) 

(602) 

(85) 


(169) 

polation 

459 

2411 


344 


678 


Table 1: Comparison of Data Storage in Number of Floating Point Data. 


The amount of floating point data in the input data and the fitting curves is 
listed in Table 1. Every input data point is counted as five floating point data: 
three for its position and two for its tangent direction, except when the data 
is in the x-y plane, in which case every input data point is counted as three 
floating point data. For the fitting curves, we count three floating point data 
for each control point, except when the curve is in the x-y plane, in which 
case we count two data for each control point. Each knot in the fitting curves 
is also counted as one floating point data. Table 1 also compares the current 
method with the rational fitting method in [9]. In the table, the number in 
the parentheses is the number of control points in the fitting curve. 


As mentioned earlier, for most of the tolerances given in Table 1, least-squares 
methods result in interpolation. Hence, the number of floating point data from 
interpolation is also listed in Table 1 for comparison. 

Table 2 lists the order of magnitude of the maximum deviation of the tangent 
direction of the fitting curves. The deviation is measured at the data points. 
Table 2 also contains the acceptance ratio, the inverse of the ratio of the 
number of closest point computation versus the number of points accepted 
without the closest point computation. The larger this ratio, the more closest 
point computations have been saved by the quick acceptance tests. 


8 























Half- 

Circle 

Torus- 

Plane 

Planar- 

Spiral 

Mixed- 

Spiral 

Tolerance 

10" 3 

IO" 6 

10" 2 

10~ 4 

2 • IO" 1 

10" 3 

IO” 5 

10’ 3 

10" 5 

Tangent 

Deviation 

10~ 3 

10" 6 

10" 2 

io - 3 

io - 1 

IO" 3 

10“ 6 

10" 3 

10' 6 

Acceptance 

Ratio 

.14 

.01 

.66 

.03 

1.27 

1.67 

.45 

.95 

2.42 


Table 2: Tangent Deviation and Acceptance Ratio. 

6 Discussion and Conclusion 


From Table 1 we can see that the current method provides much better data 
reduction than the previous methods. The difference in the reduction ratio be- 
comes wider when the tolerance is smaller. However, the current method only 
provides G 1 continuous curves. Although the fitting curve can be made C 1 by 
adjusting knot intervals, the resulting curve may have very bad parametriza- 
tion. With the rational fitting method in [8, 9], the curve can be made C 1 by 
adjusting weights of the curve while maintaining a reasonable parametric flow. 
Even though interpolation methods produce curves with the most data, they 
produce C 2 continuous curves, which may be important in some applications. 

Table 2 shows that the tangent direction of the fitting curves has roughly the 
same error as the position of the curves, even though we do not explicitly test 
the tangent direction against the given tolerance. From all the test cases, we 
see the effectiveness of the quick acceptance test is mixed. It seems to be less 
effective when the data is planar. 

Although the method presented here demonstrates great improvement in data 
reduction compared with existing methods, it still produces curves that require 
large amount of storage, especially when the required accuracy is very high. 
For example, in the Mixed-Spiral case, when tolerance = 10~ 5 , 85 control 
points are in the fitting curve for the 101 data points. Also, the storage for 
the knots becomes large, 89 knots in this case. 
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Figure 1: The given data points and tangent directions for an interpolation 
cubic. 
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Figure 2: A nonrational fit of the data set Half-Circle with 10 3 tolerance. 



Figure 3: A nonrational fit of the data set Torus-Plane with 10 2 tolerance. 
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